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ABSTRACT 

We have studied numerically the evolution of protostellar disks around intermediate and upper 
mass T Tauri stars (0.25 Mq < M* < 3.0 Mq) that have formed self-consistently from the collapse of 
molecular cloud cores. In the T Tauri phase, disks settle into a self-regulated state, with low-amplitude 
nonaxisymmetric density perturbations persisting for at least several million years. Our main finding 
is that the global effect of gravitational torques due to these perturbations is to produce disk accretion 
rates that are of the correct magnitude to explain observed accretion onto T Tauri stars. Our models 
yield a correlation between accretion rate M and stellar mass M* that has a best fit M oc , in 
good agreement with recent observations. We also predict a near-linear correlation between the disk 
accretion rate and the disk mass. 

Subject headings: accretion, accretion disks — hydrodynamics — instabilities — ISM: clouds — stars: 
formation 
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1. INTRODUCTION 

Gaseous circumstellar disks are known to extend any- 
where from several to perhaps hundreds of AU around 
T Tauri stars (TTS), due to the detection of millime- 
ter a nd submillimeter emis s ion from the associated dus t 
(e.g., iBeckwith et~aT1ll990t I Andrews fe Williams! 120051 ) . 
It has recently been recognized that such disks also ex- 
ist around brown dwarfs (BD), with masses as small 
as 0.01 Mq, but with c omparable d i sk-to-star mass ra- 
tios as found for TTS (jKlein et all 120031 : IScholz et all 
l2006f ) . Spectroscopic observations demonstrate that BD 
share similar accretion properties with TTS. In partic- 
ular, both have a large scatter in the mass accretion 
rate (2 — 3 orders of magnitud e) at a given stellar mass 
(jScholz fe Javaward hana 2006|) and, when data for BD 
and TTS are taken together, the accretion rates M show 
a strong direct dependence on stellar masses M*, with 
an approximate scaling M oc M} (e.g.. |Muzerolle et al.l 



an ap proximate scaling m oc iviz (e.g., [M uzcroiic i 
I20M 120051 iMohantv et al.ll27)MlNatta et a.l.ll200ff l 



The origin of t he above relation is uncertain. 
iPadoan et alj (|2005l ) have argued that this relation is 
a consequence of Bondi-Hoyle accretion from the large- 
scale gas distribution in the parent cloud. Such a sce- 
nario neglects the importance of disk physics in deter- 
mining the accretion rate, and also fails to explain the 
similarity of accretion rates onto TT S located within 
mole cular clouds and H II regions (jHartmann et al.l 
2006). An alternative idea is that the accretion onto 
the s tellar surface is contro l led by vis cous disk evolu- 
tion. [ Alexander fc Armitagd (|2006l ) and lHartmann et al.l 
(2006) have invoked different initial input parameters or 
scaling re lations into the standar d viscous disk evolution 
models of lHartmann et all 0.998) in order to explain the 
approximate M cx M| scaling. These models treat disks 
as isolated smooth axisymmetric structures that evolve 
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due to an unspecified source of turbulent viscosity that 
has ad-hoc spatial dependence as w ell as temporal inde- 
pendence. iDullemond et all (|2006f) have taken a more 
fundamental approach of linking the disk evolution to 
the properties of the collapsing core from which it forms, 
but still relies upon the ad-hoc a-viscosity prescription 
to model the angular momentum transport within the 
disk. 

In this Letter, we take the basic approach of study- 
ing the mass accretion rate within disks that have 
formed self-consistently due to the collapse of cloud cores. 
The disks are actually profoundly nonaxisymmetric 
when formed in this manner, and experience significant 
envelope-induced gravitational instab ility and accretion 
bursts duri ng their early evolution (jVorobvov fe Basul 
120051 120061 . hereafter VB05,VB06, respectively). Here, 
we explore the late evolution of such self-consistently 
formed disks, when their evolution is governed by 
low amplitude nonaxisy mmetric density perturbations 
(jVorobvov fc Ba su 2007, hereafter VB07). The gravita- 
tional torques produced by these perturbations are suffi- 
cient to explain the observed magnitudes and scatter of 
accretion rates of TTS, and even approximately fit the 
observed M — M* relation. 

2. MODEL DESCRIPTION 

We use the thin-disk approximation to compute the 
evolution of rotating, gravitationally bound cloud cores. 
This allows efficient calculation of the long-term evo- 
lution of a large number of models. The derivation of 
the relevant equations, details of the numerical code and 
tests are given in VB06. The basic equations of mass and 
momentum transport are 



dt 
, dv p 
' dt 
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where £ is the mass surface density, V is the vertically in- 
tegrated gas pressure, v p is the velocity in the disk plane, 
g p is the gravitational acceleration in the disk plane, and 
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V p is the gradient along the planar coordinates of the 
disk. Equations (p} and ((2]) are closed with a barotropic 
equation that makes a transition from isothermal to adia- 
batic evolution at E = S cr = 36.2 g cm~ 2 . This approach 
bypasses the detailed cooling and heating mechanisms, 
but provides a good fit to the density-temperature re- 
lation in collapsing cloud cores (see VB06). In the late 
phase of disk evolution, the stellar irradiation is expected 
to be a significant source of heating and is neglected in 
our simulation. Much of the region of enhanced tempera- 
ture would fall within our "sink cell" (see below), but our 
calculated temperatures at ~ 10 AU are also somewhat 
lower than found i n models of stellar irradia tion onto 
flared passive disks (|Chiang &: Goldreicblll997l ). Further 
details of such a comparison are given by VB07. 

Equations (fTJ) and ^ are solved in polar coordinates 
(r, (f>) on a numerical grid with 128 x 128 points. The 
radial points are logarithmically spaced. The innermost 
grid point is located at r — 5 AU, and the size of the 
first adjacent cell is 0.3 AU. We introduce a "sink cell" 
at r < 5 AU, which represents the central protostar plus 
some circumstellar disk material, and impose a free in- 
flow inner boundary condition. The outer boundary is 
such that the cloud has a constant mass and volume. 

The gas has a mean molecular mass 2.33 mn and is 
initially isothermal with temperature T (= 10 K) and 
isothermal sound speed c s . The initial distributions of 
E and angular velocity Q are those characteristic of a 
collapsing a xisymmetric magnetically supercritical core 
(|Basulll997l l: 

- (3) 
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n = 2n 
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The asymptotic r 1 power-law dependence of these 
quantities is a robust resul t of ro t ating cloud collaps e 
simulations (jNorman et al.l 119801 : iNarita et all 11984( 1. 
These profiles have the property that the specific angular 
momentum j = fir 2 is a linear function of the enclosed 
mass ro. Th e scale length r = /cc 2 /(GE ), where k ~ 1 
(|Basul [l997ll . For the models in this paper, we adopt 
k = V2/tt. These initial profiles are characterized by the 
important dimensionless free parameter 7 = V^r^/c 2 . 
The asymptotic (r > ro) ratio of centrifugal to gravita- 
tional acceleration has mag nitude ^7 (see lBasulll997ll 
and the centrifugal radius of a mass shell initially located 
at radius r is estimated to be r c f = j 2 /(Gm) — y/2jr. 
Since the enclosed mass m is a linear function of r at 
large radii, this also means that r c f oc to. 

We present results from three sets of models in this Let- 
ter, each with a different value of 7. The standard model 
has 7 = 71 = 1.2 x 10 -3 based on typical values c s = 0.19 
km s _1 , E = 0.12 g cm -2 , and f^o = 1.0 km s _1 pc" 1 . 
The outer radius is taken to be r ou t = 0.04 pc, and the 
total cloud mass is 0.8 Mq. Other models with 7 = 71 
but different mass are generated by varying ro and fio 
so that their product is constant. All clouds are char- 
acterized by the same ratio r ou t/^o ~ 6.0. To generate 
the second set of models, with 7 = 72 = 2.3 x 10~ 3 , we 
set J7q = 1-4 km s _1 pc -1 and all other quantities the 



same as in the standard model with 7 = 71. Models 
of varying mass are then generated in the same man- 
ner as for the 71 models. The third set of models, with 
7 = 73 = 3.4 x 10 -3 , are also obtained in this way, by 
first using f2o = 1-7 km s _1 pc -1 . Overall, there are 6 
models with 7 = 71, 14 models with 7 = 72 ~ 271, and 
12 with 7 = 73 ~ 371. The range of initial cloud masses 
amongst our models is 0.25 Mq — 3.0 Mq. 

The numerical simulations start in the prestellar phase 
and continue into the late accretion phase, long after the 
formation of a protostar and a protostellar disk. The 
disk evolution is followed for approximately 3 Myr after 
the formation of the protostar. In the early phase, when 
the infall of matter from the surrounding envelope is sub- 
stantial, mass is transported inward by the gravitational 
torques from spiral arms that are a manifestation of 
the envelope-induced gravitational instability in the disk 
(VB05; VB06). In the late phase, when the gas reservoir 
of the envelope is depleted, the distinct spiral structure 
is replaced by ongoing irregular nonaxisymmetric density 
perturbations in the disk. These perturbations are con- 
fined to the disk, of size ~ 100 AU, and are not affected 
by the computational boundary at ~ 5, 000—10, 000 AU. 
Instead, their longevity is enhanced by swing amplifica- 
tion at the disk's sharp physical boundary (see VB07). 
We find that the net effect from these density pertur- 
bations is a residual non-zero gravitational torque. In 
particular, the net gravitational torque in the inner disk 
tends to be negative during first several million years 
of the evolution, while the outer disk has a net posi- 
tive gravitational torque. The inward radial transport of 
matter due to the negative torque, which is produced self- 
consistently in our numerical hydrodynamic modeling, is 
the essence of the accretion mechanism in our model. 

Our model accretion rates are consistent with typi- 
cal accretion rates for intermediate and upper-mass TTS 
(0.25 Mq < M* < 3.0 M Q ). We do not consider objects 
with masses below 0.25 Mq, because the numerical noise 
generated by the inner boundary becomes comparable to 
the amplitude of density perturbations in compact disks 
around extremely low-mass objects (< 0.1 Mq). 

3. ACCRETION RATES 

Figure[T]shows confirmed detections for the mass accre- 
tion rate M and stellar mass M* for TTS and BD of age 
0.5 — 3 Myr from two recent observational compilations. 
The open diamonds represent measure ments, mostly in 
Tauru s, that have been compiled by iMuzerolle et al.l 
((20051 and references therein); the open squares repre- 
sent detections in p Oph obtained bv lNatta et al.l (2006) 
and crosses represent their upper limits to nondetec- 
tio ns. We have exclude d objects in the compilation 
of IMuzerolle et al l (|2005( 1 that were later observed by 
iNatta et al.l (|2006H . The least-squares best fit to the ob- 
servational data, both TTS and BD (upper limits ex- 
cluded), has an exponent 2.0 ±0.1 and is shown in Fig- 
ure[T]by a black line. This value is often quoted in the lit- 
erature (see § 1). However, we believe that taking a best 
fit over the whole mass range of BD and TTS may be mis- 
leading. Indeed, if we consider separately the intermedi- 
ate and upper-mass TTS (0.25 Mq < M* < 3.0 Mq) and 
lower-mass TTS plus BD (0.02 M < M, < 0.25 Mq), 
then the least-squares best fits to each data sample are 
distinct. In particular, for lower-mass TTS and BD we 
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Fig. 1. — Mass accretion rates versus stellar masses. The open diamonds represent measurements of TTS and BD from Muzcrollc ct al. 
(120051. and references therein); the open squares and crosses represent the confirmed detections and upper limits, respectively, compiled 
by INatta et~aTl &200(i) . The filled triangles, filled circles, and filled squares show the time-averaged quantities for the models with 7 = 
71,72, and 73, respectively. The bars represent mean positive/negative deviations from the time-averaged accretion rates in each model. 
The blue line (oc ) is the least-squares fit to the model data. The black line (oc Mj'°) is the least-squares fit to the observed confirmed 
detections, both TTS and BD. The left and right red lines are the least-squares fits to objects with masses in the range (0.02 — 0.2) Mq 
and (0.2 — 3.0) Mq, respectively. 



obtain an exponent 2.3 ± 0.6 (left red line), whereas for 
the intermediate and upper-mass TTS we find a much 
smaller exponent 1.3 ± 0.3 (right red line). 

The above data hints that different mechanisms may be 
responsible for accretion as one moves along the sequence 
of stellar masses. However, the observational method of 
determining M also typically varies across the mass se- 
quence, with significant but differing uncertainties. For 
TTS, the primary determinant of M at the stellar surfac e 
is UV excess and veiling (see, e.g. lMuzerolle et aTi r2003'). 
while for BD t he primary means is th e fitting of emission 
line profiles. iMuzerolle et alj (|2003l ) point out that the 
latter is considered less accurate than UV excess mea- 
surements; it suffers from uncertainties regarding optical 
depth, rotation, and inclina tion effects (see also discus- 
sion in lMohantv et all 12005). Quantitative estimates of 
error bars for either technique are not available in the lit- 
erature. However, the agreement between the two meth- 
ods is within a factor ^3 — 5 (|Muzerolle et al.l l2003f) 
where comparable, and this is much smaller than the 
spread of observed M for a given central object mass, 
implying that the spread of values is primarily due to 
physical effects of disk age and initial conditi ons. We 
also note that the data set of lNatta et alj (|2006f ) is based 
entirely on emission line estimates of accretion onto BD 
and TTS. 

Figure Q] also shows the time-averaged mass accretion 
rates (M) and time- averaged stellar masses (M*) for our 
models with 7 = 71,72, and 73, respectively. The mass 



accretion rate M(t) = —2nrv r Y,, where v r is the inflow 
velocity of gas through the sink cell and r = 5 AU. For 
most models, the time average is taken between 0.5 Myr 
and 3.0 Myr after the formation of the protostar. How- 
ever, protostars with masses above 2.0 Mq may enter 
the T Tauri phase (class II) when they are older than 
0.5 Myr. To exclude a possible input from class 0/class I 
sources in such cases, we start the time average only when 
the envelope mass has dropped below 10% of the initial 
cloud mass. The best fit to our model data points is 



(M) = 10~ 7 ' 7 (M*) 



(5) 



and is represented in Figure [T] with a blue line. The 
least-squares method generates an uncertainty ±0.1 in 
the above exponent. The somewhat shallower best-fit 
slope to the observational data (1.3±0.3, right red line) is 
likely because the short-lived FU Ori bursts are not sam- 
pled observationally given the small number of objects. 
On the other hand, our numerical models produce FU- 
Ori-like mass accretion bursts (see Fig. [3]), which effec- 
tively increase (M) and steepen the model best-fit slope. 

It should be noted that our model accretion rates are 
derived at 5 AU (the sink cell), while the observed accre- 
tion rates are measured near the stellar surface. However, 
our numerical simulations indicate that accretion rates, 
averaged over many orbital periods, vary little with ra- 
dius in the inner several tens of AU. Some physical pro- 
cesses operating in the inner several AU and unaccounted 
in our numerical modeling, will allow accretion onto the 
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Fig. 2. — Time-averaged mass accretion rate (M) versus time- 
averaged disk mass (Afj) for all models. The solid line shows the 
least-squares fit to the data points. Symbols have the same mean- 
ing as in Fig. [TJ 

stellar surface. The inner accretion rate is expected to 
match our calculated value when time-averaged, even 
though it may have significant short term variability. 

To quantify the characteristic range of accretion rates 
obtained in our models, Figure Q] also shows, for each 
model, the upper and lower bounds on the mass accre- 
tion rate M. These values are obtained by smoothing 
the accretion rates over 10 yr periods (to reduce noise) 
and typically correspond to the accretion rates at the be- 
ginning of the T Tauri phase (upper bound) and at the 
terminal point of the simulations (lower bound) . It is ev- 
ident that our models can account for the observed range 
of accretion rates for TTS with masses above 1.0 Mq. On 
the other hand, the observed range of accretion rates for 
the intermediate-mass TTS (0.25 M Q < M* < 1.0 Mq) 
is greater than that implied by the temporal evolution of 
our models. We believe that this can be accommodated 
ultimately by a broader range of initial cloud configura- 
tions than we have studied here due to numerical limi- 
tations. However, it is remarkable that our models do 
cover the middle portion of the observed M — M* phase 
space. 

Figure [5] shows the relation between the time-averaged 
mass accretion rate (M) and the time-averaged disk mass 
(Md) , to which the least-squares best fit is 

(M) = 10- 7 (Md) 1 ' 1 . (6) 

The disk masses are that of matter with surface den- 
sity above a value 0.1 g cm -2 that typically charac- 
terizes the disk-to-envelope transition (VB07). In our 
view, equation ^ is the most physically meaningful cor- 
relation arising from our simulations of self-consistent 
disk formation and evolution due to global gravitational 
torques. However, there is also a mild preference for 
more massive disks around more massive stars. A plot 
of £ = (M d )/(M*) versus (M*) has a best-fit £ oc 
(M*)°- 3±01 , with values ranging from ~ 5% at the low 
mass end to ~ 35% at the high mass end. This serves to 
steepen the correlation in Figure Q] so that the best fit is 
(M) oc (M*) 1 - 7 . 

Finally, Figure[3]shows the typical mass accretion rates 
obtained in our models as a function of time. The top 
and bottom lines correspond to models with initial cloud 
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Fig. 3. — Mass accretion rate M versus time for two models, 
both with 7 = 72 and initial cloud masses 2.75 Mq (top) and 
0.35 Mq (bottom). The solid lines show the least-squares best fits 
to the data points. 

masses 2.75 Mq and 0.35 Mq, respectively, both with 
7 = 72. The corresponding time- average protostellar 
masses are 0.63 Mq and 0.026 Mq, respectively. The 
more massive disks clearly drive greater mass accretion 
rates. Furthermore, more massive disks show a wider 
range of mass accretion rates in the T Tauri phase. The 
least-squares best fit to the data in Figure [3] yields an 
exponent —3.5 for the (Md) = 0.63 Mq disk and —1.2 
for the (Md) = 0.026 Mq disk. Massive disks have grav- 
itational torques of greater magnitude, which result in 
greater associated accretion rates. The more massive 
disk is also violently gravitationally unstable in its early 
evolution and is characterized by FU-Ori-like accretion 
bursts (VB05; VB06). We believe that the steeper de- 
cline of accretion rate of the (M d ) = 0.63 Mq disk is 
caused by the greater effect of disk self-gravity and pos- 
sibly the influence of the more massive stellar object at 
the center. 

4. CONCLUSIONS 

In this Letter, we have presented models of the late evo- 
lution of self-consistently formed protostellar disks that 
can explain the observed correlation of disk accretion 
rates with stellar mass for TTS. The formation of the 
disk and its interaction with the surrounding envelope 
lead to the development of strong spiral structure during 
the early evolution of the disk (VB05; VB06). Even after 
the former has largely disappeared, low-amplitude non- 
axisymmetric density perturbations are sustained in the 
disk for several Myr. The gravitational torques due to 
these perturbations are sufficient to drive accretion at the 
rates commonly inferred around TTS. We find that an 
important property of gravitational torques is that (M) 
has an essentially linear dependence on (M d ). The aver- 
age disk-to-star mass ratio £ is in the range ~ 5% — 35% 
for models of various masses, given our adopted range 
of values of initial cloud angular momentum. However, 
there is a mild trend toward greater values of £ for models 
with greater masses. The net result is a best-fit corre- 
lation (M) oc (A/*) 1-7 , in good agreement with the ob- 
served M — M* relation. 

Since the of values of £ in our models are ~ 10 times 
greater than estimates made from dust emission (e.g., 
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lAndrews fe Williamsl 12003 iScholz et al.ll200l . wc an- 
ticipate two reasons that this discrepancy may be re- 
duced in the future. Observationally, the inferred gas 
disk masses may be systematically un derestimated us- 
ing c urrent techniques (see discussion in lHartmann et all 
2006). Theoretically, wc need to include additional angu- 
lar momentum transport mechanisms such as magnetic 
braking and the magnetorotational instability. We note 
that only a modest increase in the overall average accre- 
tion rate in our current models is required to lead to a 
dramatic decrease in £. This is because a small relative 
increase in the stellar mass can significantly reduce the 
much smaller disk mass. 
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